function generateMfuncResidualEquationBondPrices_file(allModelParams,x,y,xp,yp,nx1,ne,h11,heteroShock,BondEq,positionFile,nameOfFunction)
% This function generates an M function which is the objective function
% when solving for bond prices by projection

% We start deleting the old version of the file - if it exists
saveFile = [positionFile,'\',nameOfFunction];
if exist(saveFile,'file') > 0
    delete(saveFile)
end

text = ['function [res,resEuler,p2_cup,logP2_cu] = ',nameOfFunction(1:end-2),'(thetaVec,setup)'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%% Getting some information from setup ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'scaleX_in_g = setup.scaleX_in_g;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'out         = setup.out;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'x_cu        = setup.xGrid;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'numX        = size(x_cu,1);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nNodes      = setup.n_nodes;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'wNodes      = transpose(setup.weight_nodes);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'epsNodes    = setup.epsi_nodes;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'params      = setup.params;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'appOrder    = setup.appOrder;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nx          = setup.nx;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'ny          = setup.ny;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nx1         = setup.nx1;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'gSS         = setup.gSS;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'hSS         = setup.hSS;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% For the new bond price - the unknowns';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'thetaProj   = reshape(thetaVec,1,setup.nDim2Theta);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'g0          = thetaProj(1,1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'gx          = thetaProj(1,2:1+nx); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    gxxV    = thetaProj(1,nx+2:1+nx+nx*(nx+1)/2); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 2 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    gxxxV   = thetaProj(1,1+nx+nx*(nx+1)/2+1:1+nx+nx*(nx+1)/2+nx*(nx+1)*(nx+2)/6); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%Unfold params ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(allModelParams)
    text = [allModelParams{i}, '= params.',allModelParams{i},';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end 
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%% Unfold the states and controls for the core of the model - all in levels ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the states at time t ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
     text = [char(x(i)),' = out.',char(x(i)),';'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the controls at time t ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    text = [char(y(i)),' = out.',char(y(i)),';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the states at time t+1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(xp)
     text = [char(xp(i)),' = out.',char(xp(i)),';'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the controls at time t+1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(yp)
     text = [char(yp(i)),' = out.',char(yp(i)),';'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Bond price ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'p1_cup = out.p1_cup;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Computing the new bond price ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
    text = ['x',num2str(i),'_cu ',' = x_cu(:,',num2str(i),')/scaleX_in_g(',num2str(i),',1);'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
name = 'p2_cu = g0(1,1)';
for j=1:length(x)
    name = [name,['+gx(1',',',num2str(j),').*x',num2str(j),'_cu']];
end
text = [name,';'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 1';
dlmwrite(saveFile,text,'-append','delimiter','');
idx = 0;
name = 'p2_cu = p2_cu';
for alfa1=1:length(x)
    for alfa2=alfa1:length(x)
        idx = idx + 1;
        if alfa1 == alfa2
            name = [name,['+gxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu']];
        else
            name = [name,['+2.*gxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu']];
        end
    end
end
text = ['    ',name,';'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 2';
dlmwrite(saveFile,text,'-append','delimiter','');
name = 'p2_cu = p2_cu';
idx = 0;
for alfa1=1:length(x)
    for alfa2=alfa1:length(x)
        for alfa3=alfa2:length(x)
            idx = idx + 1;
            if alfa1 == alfa2 && alfa2 == alfa3
                name = [name,['+gxxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
            elseif alfa1 == alfa2 && alfa2 ~= alfa3
                name = [name,['+3.*gxxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
            elseif alfa1 ~= alfa2 && alfa2 == alfa3
                name = [name,['+3.*gxxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
            elseif alfa1 ~= alfa2 && alfa2 ~= alfa3
                name = [name,['+6.*gxxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
            else
                [alfa1,alfa2,alfa3];
            end
        end
    end
end
text = ['    ',name,';'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% The bond price next period ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if setup.computeP2_cup == 1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    xVec_cup = zeros(numX,nx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    mx = setup.mx;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    if mx > 0 ';
dlmwrite(saveFile,text,'-append','delimiter','');
nx11 = length(h11);
for i=1:nx11
    TEMP = char(h11(i,:));
    TEMP = strrep(TEMP,'*','.*');
    TEMP = strrep(TEMP,'/','./');
    TEMP = strrep(TEMP,'^','.^');
    text = ['    ',char(xp(1,i)),' = ',TEMP,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = '    end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    if setup.myx > 0 ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:nx1-nx11
    text = ['        xVec_cup(:,mx+',num2str(i),') = ',char(y(1,i)),'(:,1)-gSS(',num2str(i),',1);'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = '    end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    diaghx = diag(setup.hx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    xVec_cup(:,nx1+1:nx)  = x_cu(:,nx1+1:nx).*repmat(transpose(diaghx(nx1+1:nx)),numX,1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    %Adding shocks to the states at t + 1';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    diagEta = diag(setup.eta(nx1+1:nx,:));';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:nx1
    name = ['x',num2str(i),'_cup  = repmat(xVec_cup(:,',num2str(i),'),1,nNodes);'];
    text = ['    ',name];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
for i=1:ne
    name = ['x',num2str(nx1+i),'_cup  = repmat(xVec_cup(:,',num2str(nx1+i),'),1,nNodes) +',...
         '2./(exp(-',char(heteroShock(i,1)),'.*x_cu(:,',num2str(nx1+i),'))+1).*diagEta(',...
         num2str(i),',1).*repmat(epsNodes(',num2str(i),',:),numX,1);'];
    text = ['    ',name];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    % Scaling the states at time t+1 by the stdX';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
    name = ['x',num2str(i),'_cup  = x',num2str(i),'_cup/scaleX_in_g(',num2str(i),',1);'];
    text = ['    ',name];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
name = 'p2_cup = g0(1,1)';
for j=1:length(x)
    name = [name,['+gx(1',',',num2str(j),').*x',num2str(j),'_cup']];
end
text = ['    ',name,';'];
dlmwrite(saveFile,text,'-append','delimiter','');

text = '    if appOrder > 1';
dlmwrite(saveFile,text,'-append','delimiter','');
name = 'p2_cup = p2_cup';
idx = 0;
for alfa1=1:length(x)
    for alfa2=alfa1:length(x)
        idx = idx + 1;
        if alfa1 == alfa2
            name = [name,['+gxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup']];
        else
            name = [name,['+2.*gxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup']];
        end
    end
end
text = ['        ',name,';'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    if appOrder > 2';
dlmwrite(saveFile,text,'-append','delimiter','');
idx = 0;
name = 'p2_cup = p2_cup';
for alfa1=1:length(x)
    for alfa2=alfa1:length(x)
        for alfa3=alfa2:length(x)
            idx = idx + 1;
            if alfa1 == alfa2 && alfa2 == alfa3
                name = [name,['+gxxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
            elseif alfa1 == alfa2 && alfa2 ~= alfa3
                name = [name,['+3.*gxxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
            elseif alfa1 ~= alfa2 && alfa2 == alfa3
                name = [name,['+3.*gxxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
            elseif alfa1 ~= alfa2 && alfa2 ~= alfa3
                name = [name,['+6.*gxxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
            else
                [alfa1,alfa2,alfa3];
            end
        end
    end
end
text = ['        ',name,';'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'else';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    p2_cup = NaN;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% The Euler errors ';
dlmwrite(saveFile,text,'-append','delimiter','');
TEMP = BondEq;
TEMP = strrep(TEMP,'*','.*');
TEMP = strrep(TEMP,'/','./');
TEMP = strrep(TEMP,'^','.^');
text = ['E_sdf_p1_cup = ',TEMP,';'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'f1 = E_sdf_p1_cup./exp(p2_cu) - 1;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Numerical integration ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'resEuler = f1*wNodes;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'logP2_cu = log(E_sdf_p1_cup*wNodes);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'resEulerWeight = resEuler.*(1./(1+sum(x_cu.^2,2)));';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'res = resEulerWeight(:);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
end